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We study the problem of pattern selection in an array of parametrically-driven nonlinear res- 
onators with application to microelectromechanical and nanoelectromechanical systems (MEMS & 
NEMS), using an amplitude equation recently derived by Bromberg, Cross, and Lifshitz [Phys. Rev. 
E 73, 016214 (2006)]. We describe the transitions between standing- wave patterns of different wave 
numbers as the drive amplitude is varied either quasistatically, abruptly, or as a linear ramp in 
time. We find novel hysteretic effects, which are confirmed by numerical integration of the original 
equations of motion of the interacting nonlinear resonators. 



PACS numbers: 45.70.Qj, 62.25.-g, 85.85.-fj, 05.45.-a 



I. INTRODUCTION 

Interest in the nonlinear dynamics of microelectrome- 
chanical and nanoelectromechanical systems (MEMS & 
NEMS) has grown rapidly over the last few years, driven 
by a combination of practical needs as well as fundamen- 
tal questions [l[ . Lithographic fabrication technology al- 
lows the construction of large arrays of MEMS & NEMS 
devices (as many as 2800 to date [2]), coupled by electric, 
magnetic, or elastic forces. In addition, nonlinear behav- 
ior is readily observed in these devices at relatively small 
amplitudes of motion fl, i, i, i, 0, S S [M M MM,M - 
Limitations in the fabrication technology mean that indi- 
vidual devices will usually have slightly different resonant 
frequencies, and nonlinear collective effects, such as syn- 
chronization (all devices oscillating in phase) [15, 16] and 
pattern formation [13, [H, Elj HSI (coherent response with 
a more complex spatial structure), have been proposed 
as ways of achieving useful coherent responses. Conse- 
quently, for many technological applications, there exists 
a practical need to understand the collective nonlinear 
behavior of MEMS & NEMS devices. 

At the same time, the advances in the fabrication, 
transduction, and detection of MEMS & NEMS res- 
onators opens up an exciting new experimental win- 
dow into the study of fundamental questions in collec- 
tive nonlinear dynamics. Typical nonlinear MEMS & 
NEMS resonators are characterized by extrernely high 
frequencies — recently going beyond 1 GHz [2l|, — and 
relatively weak dissipation, with quality factors in the 
range of 10^ — 10^. For such devices, transients die out 
rapidly, so that it is easy to attain the long-time asymp- 
totic states, be they steady, periodic, or chaotic, and to 
acquire sufficient data to characterize these states well. 
From the theoretical point of view, the systems have the 
advantage that the basic physics of the individual ele- 
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ments is simple, and the parameters can be measured or 
calculated, so that the equations of motion describing the 
system can be established with confidence. This, and the 
fact that weak dissipation can be treated as a small per- 
turbation, provide a great advantage for quantitative the- 
oretical study. Moreover, the ability to fabricate arrays 
of thousands of coupled resonators opens new possibili- 
ties in the study of nonlinear dynamics of intermediate 
numbers of degrees of freedom — much larger than one 
can study in macroscopic or table-top experiments, yet 
much smaller than one studies when considering nonlin- 
ear aspects of phonon dynamics in a crystal. 

Our current studies are motivated by the experimen- 
tal work of Buks and Roukes [13] , who fabricated an ar- 
ray of nonlinear micromechanical doubly-clamped gold 
beams, and excited them parametrically by modulating 
the strength of an externally-controlled electrostatic cou- 
pling between neighboring beams. The Buks and Roukes 
experiment was modeled by Lifshitz and Cross [l^ us- 
ing a set of coupled nonlinear equations of motion. They 
used secular perturbation theory to convert these equa- 
tions of motion into a set of coupled nonlinear algebraic 
equations for the normal mode amplitudes of the system, 
enabling them to obtain exact results for small arrays, 
but only a qualitative understanding of the dynamics of 
large arrays. In order to obtain analytical results for 
large arrays, Bromberg, Cross, and Lifshitz [13, hence- 
forth BCL] studied the same system of equations, ap- 
proaching it from the continuous limit of infinitely-many 
degrees of freedom, and obtaining a description of the 
slow spatiotemporal dynamics of the array of resonators 
in terms of an amplitude equation. BCL showed that 
this amplitude equation could predict the initial mode 
that develops at the onset of parametric oscillations as 
the driving amplitude is gradually increased from zero, 
as well as a sequence of subsequent transitions to other 
single-mode oscillations. 

The combination of many degrees of freedom and non- 
linearity in the equations of motion typically leads to 
a large multiplicity of physically realizable solutions for 
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fixed system parameters. This is illustrated for the par- 
ticular case of two and three parametrically driven oscil- 
lators by the explicit results of Lifshitz and Cross p^ . 
The richness of possible solutions leads to opportunities 
for diverse functionality of the system, in nature or tech- 
nology. On the other hand we need to be able to pre- 
dict which out of the possible solutions will be seen for 
a given experimental protocol, or design particular pro- 
tocols such that the desired solution is the one that is 
formed. This is the general question of pattern selec- 
tion (23j . A common experimental protocol is to vary 
one or more system control parameters, usually either 
slowly compared with the intrinsic time scales of the dy- 
namics, or in an abrupt step. A particular solution will 
usually survive (evolving adiabatically in the former case 
of slow parameter variation) until it becomes unstable to 
small perturbations, and a sequence of patterns can be 
predicted by analyzing these instabilities. 



In this paper we investigate the sequence of single 
mode standing wave patterns to be expected in paramet- 
rically driven oscillator arrays, in cases where many such 
modes are simultaneously stable, when the strength of 
the driving is varied. Although the quantitative analysis 
could be done directly from the basic oscillator equations 
of motion, it is advantageous to formulate the analysis in 
terms of the BCL amplitude equation. This allows us 
to display the range of stable patterns on a reduced plot 
involving just two dimensionless variables (a scaled mea- 
sure of the driving strength, and a scaled mode wave 
number), so that it is easy to deduce the general qualita- 
tive behavior on varying parameters. The specific quan- 
titative behavior for a physical system is also easy to ob- 
tain by evaluating the corresponding scaled quantities. A 
change of pattern occurs when parameters vary so that 
the mode moves outside of the region of stable patterns 
on this plot, and the new pattern is predicted by ana- 
lyzing the result of the instability using the BCL am- 
plitude equation. This type of approach has been used 
in other pattern forming systems 12%- A novel feature 
of the present system is that the difference in the in- 
stabilities encountered on increasing and decreasing the 
(scaled) driving strength leads to the prediction of quite 
different sized mode jumps for the up and down sweeps. 



The outline of the paper is as follows. In Sec. [IT] we 
review the derivation of the amplitude equation of BCL, 
and in Sec. IIIII use this equation to discuss the stability 
of single-mode oscillating patterns. We then study the 
sequence of patterns observed for a variety of time depen- 
dent sweeps of the driving strength: quasistatic variation 
in Sec. lIVI abrupt step jumps in Sec.|Vl and a control pa- 
rameter ramp varying linearly in time in Sec. lVIl Finally, 
we conclude with some remarks connecting our results to 
those of Buks and Roukes [13] who swept the frequency 
rather than the driving strength. 



II. BCL AMPLITUDE EQUATION 

Lifshitz and Cross [l^ modeled the array of cou- 
pled nonlinear resonators that was studied by Buks and 
Roukes [13] using the equations of motion 

■Un+Un+ ul - ^e{Un+l - 2u„ + Un^l) 

+ i [A^ -I- ehcos{2Ljpt)] (u„+i - 2u„ -|- Un-i) 

- ^'^^^^[("n+l - U„)^(m„+1 - Un) 

" (u„ - u„_i)^(m„ - U„_l)] = 0, (1) 

where Un{t) describes the deviation of the n*^ resonator 
from its equilibrium, with n = 1 . . . N, and fixed bound- 
ary conditions uq = uat+i = 0. Detailed arguments for 
the choice of terms introduced into the equations of mo- 
tion are discussed in Ref. [l^. The terms include an 
elastic restoring force with both linear and cubic con- 
tributions (whose coefficients are both scaled to 1), a dc 
electrostatic nearest-neighbor coupling term with a small 
ac component responsible for the parametric excitation 
(with coefficients and eh respectively), and linear as 
well as cubic nonlinear dissipation terms. The dissipation 
in the system is assumed to be weak, which is used to de- 
fine two small expansion parameters e <C 1 and (5 <SC 1 
by setting the linear damping rate to e and the nonlinear 
damping coefficient to S^^^, with a square root for later 
convenience. The driving amplitude is then expressed 
as eh, with h of order one, in anticipation of the fact 
that parametric oscillations at half the driving frequency 
require a driving amplitude which is of the same order 
as the linear damping rate (25j . Both dissipation terms 
are taken to be of a nearest neighbor form, motivated by 
the experimental indication that most of the dissipation 
comes from the electrostatic interaction between neigh- 
boring beams. 

In order to treat the system of equations ^ analyt- 
ically, BCL introduced a continuous displacement field 
u(x,t), and slow spatial and temporal scales, X = ex 
and T = et. They tried a solution in terms of a pair of 
counter-propagating plane waves, oscillating at half the 
drive frequency, 

u{x,t) = e^^^[{A+{X,T)e-''''''= + A*_{X,T)e"''-'')e"^'-* 
+ c.c] + e3/2u(i) (x, t,X,T) + ..., (2) 

where the asterisk and c.c. denote complex conjugation, 
and Qp and ujp are related through the dispersion relation 

c.^ = l-2A^sin2(|). (3) 

By substituting this ansatz Q into the equations of mo- 
tion IT]) and applying a solvability condition on the terms 
of order e^/^, BCL obtained a pair of coupled amplitude 
equations for the counter-propagating wave amplitudes 
A±. A linear analysis of these equations shows that at 
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the critical drive amplitude he = a particular lin- 

ear combination of the two counter-propagating waves 
obtains a positive growth rate, forming a standing wave 
pattern, while the growth rate of the orthogonal linear 
combination remains negative. This implies that a sin- 
gle amplitude equation should suffice at onset, describing 
this standing wave pattern. 

At this point it is natural to define a reduced driving 
amplitude g with respect to the critical drive at onset 
by letting {h — hc)/hc = gS, and to introduce a second 
ansatz. 



0' ]]B{lf) + S'/' 



,(1) 



(^,T,e,T) 



5/4 / w^^\X,T,if) 



(4) 



where ^ = S^^^X and f = ST. Substitution of this ansatz 
allows one to obtain the correction of the solution at or- 
der 53/4 



yd) 
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Aujp sin^(gp/2) 
dB 



X sin (qp) — + 9i\B\'^B 
\ 9^ ) 



(5) 



after which a solvability condition applied to the terms of 
order J^/"*, and a rescaling of all the physical quantities, 
yield an equation for the scaled field B(^^^ r) of the form 



dB 

a7 



9B 



d^B 



4IBI 



,95 



B 



,dB* 



2\B?B 



(6) 



This is the BCL amplitude equation. It is governed by 
a single control parameter, the reduced drive amplitude 
g, and captures the slow dynamics of the coupled res- 
onators just above the onset of parametric oscillations. 
The reader is encouraged to consult Ref. [l^l for a more 
detailed account of the derivation of the BCL equation, 
as well as a detailed list of all the scale factors leading to 
the final form of the equation. 



III. SINGLE-MODE SOLUTIONS OF THE BCL 
AMPLITUDE EQUATION 

The simplest nontrivial solutions of the BCL amplitude 
equation are steady-state single- mode extended patterns, 
given by 



B{^,T)^bke'(^-''^\ 



(7) 



with bk and Lp both real. This solution, when substituted 
back into (O and (H]), and then into yields single- 
mode standing-wave parametric oscillations at half the 
drive frequency, whose explicit form is given in Appendix 



A. The original boundary conditions M(0,i) = u{N + 
= constrain the phase ip to be 7r/4 or 57r/4, and 
constrain the wave numbers of the spatial pattern to have 
the quantized values of q„i — mn/^N -1-1), with m — 
1,...,N. 

BCL showed that the first single-mode pattern to 
emerge as the zero-state becomes unstable is that whose 
wave number is closest to the wave number Qp that is 
determined by the drive frequency ujp through the disper- 
sion relation ^ . This determines the value of the scaled 
wave number in the single-mode solution ([7]) to be 



fco = [m- Qp 



N+1 



AQ 



N, 



(8) 



where m is the integer closest to qp{N + l)/7r, and AQn, 
whose explicit value is given in Appendix A, tends to zero 
as the size N of the array of resonators tends to infinity. 
In this paper we are interested in secondary transitions as 
the initial single-mode state of wave number fco becomes 
unstable with respect to the growth of other single-mode 
states, whose wave numbers we label as 



ko + nAQ 



N- 



(9) 



In steady state, the relation between the magnitude 
and the wave number k is found by substituting ([7]) into 
([6]) and setting the time derivative to zero, to give 



bl = (fc - 1) + ^(fc - 1)2 -t- (5 - fc2) > 0, 



(10) 



along with a negative square-root branch which is always 
unstable against small perturbations (201] . as can be ver- 
ified by the analysis below. Linearization of the BCL 
amplitude equation ([5]) shows that the zero state with 
B{£,,t) = — which is a solution of ^ for any value of 
g — is stable against the formation of single-mode pat- 
terns with wave number k as long as g < k^. The neutral 
stability curve g — is plotted as a dashed parabola in 
Fig. [TJ Furthermore, for A; < 1 the bifurcation from the 
zero state to that of single-mode oscillations is supercrit- 
ical, occurring on the neutral stability curve, while for 
/c > 1 it is subcritical, with a locus of saddle-node bifur- 
cations located along the line g = 2k—l (shown in Fig.[T] 
as a solid green line), where the square root in (fTO]) is 
exactly zero. 

The stability of a single-mode solution ([7]) of wave 
number k against an Eckhaus transition to a different 
single-mode solution of wave number k ± Q is found by 
performing a linear stability analysis of solutions of the 
form 

B{^,t) = 5fce-^^-«-f (/3+(r)e-'('=+«)« +/3:i(T)e-^('=-«)«) , 

(11) 

with <C 1. When the larger of the two eigenvalues 
describing the growth of such a perturbation, which is 
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given by [lO] 

X,AQ) = 2blik-l-bl)-Q' 

+ l[3Q\k-bl){3k~5bl) 

+ 9bUk ~ 1 - bl)'Y^\ (12) 

becomes positive the single-mode solution of wave num- 
ber k undergoes an Eckhaus instability with respect to 
different single- mode solutions of wave numbers fciQfll]. 

For an infinite number of oscillators the Eckhaus insta- 
bility forms the upper boundary of the stability balloon 
of the single-mode solutions, and also the lower bound- 
ary for k < 5/2. For fc > 5/2 the lower boundary is 
the saddle node bifurcation line. For a finite number 
of oscillators, restricting Q to be an integer multiple of 
AQn in slightly shifts the Eckhaus instability lines. 
The upper Eckhaus boundary is shifted to larger val- 
ues of g. The nature of the lower instability boundary 
now depends on the number of resonators in the array 
through AQtv, as well as on the wave number k [20|. 
For fc < 1 the lower boundary will be the Eckhaus in- 
stability curve if \k\ > AQn/2, and the neutral stabil- 
ity curve otherwise. From ([5]) and ^ we find that the 
only wave number to satisfy \k\ < AQn/2 is feg, which 
means that upon decreasing g the ko solution undergoes 
a continuous transition to the zero state. For k > 1 
the lower boundary will be the Eckhaus instability curve 
if 1 < A: < (5 - 3(AgAr/2)2)/2, and the line of saddle 
node bifurcations otherwise. For AQn > 2 (for the pa- 
rameters used throughout this paper this corresponds to 
N < 172) there is no portion of Eckhaus instability on 
the lower boundary, which is the neutral stability curve 
if fc < 1 and the saddle node bifurcation curve if fc > 1. 
These stability boundaries are shown in Fig. [1] for an infi- 
nite system and for a system of iV = 92 resonators, which 
is discussed next. 



IV. QUASISTATIC SWEEPS OF THE 
CONTROL PARAMETER 

We begin by taking a close look at the switching that 
occurs between single-mode patterns ([7]) of different wave 
numbers fc„ as the control parameter — the reduced drive 
amplitude g — is varied quasistatically. We examine a 
typical situation, which is depicted within the stability 
balloon of single-mode solutions, shown in Fig. [1] Pa- 
rameters are chosen such that the initial pattern happens 
to have a wave number ko ~ —0.81, which corresponds 
to the array of = 92 nonlinear resonators oscillating at 
its TO = 67 mode. Because fcg < 1 we expect the pattern 
to grow supercritically from the zero state as the control 
parameter is gradually increased from g ~ 0. The se- 
quence of expected secondary transitions to single-mode 
patterns of wave numbers fc„ can be understood with the 
help of the vertical and horizontal lines drawn within the 
stability balloon. As g reaches a value of about 10, the 
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FIG. 1: (Color online) Stability boundaries of the single-mode 
solution (O of the BCL amplitude equation ^ in the g vs. k 
plane. Dashed line: neutral stability curve g = . Dotted 
line: stability boundary of the single-mode solution ([7]) for a 
continuous spectrum (Q ^ 0). Solid lines: stability boundary 
of the single-mode solution for TV = 92 and the parameters 
A = 0.5, qp = 737r/101, and e = 5 = 0.01 (giving ko ~ -0.81 
and AQn ~ 3.70). Black line: the value of g for which the 
eigenvalue Ag^fc(AQ]v) turns positive. Red line: the lower 
bound for k < 1, g = . Green line: the lower bound for 
fc > 1, the locus of saddle-node bifurcations g — 2k — 1. 
Vertical and horizontal arrows mark the secondary instability 
transitions shown in Fig. [2] and discussed in Sec. IIVI 



initial fco pattern undergoes an Eckhaus instability to a 
pattern of wave number fci ~ 2.90. As this occurs in the 
solution ^ of the amplitude equation ([6]) the pattern of 
the array of nonlinear resonators U) switches from the 
67*'* mode to the 68*'' mode via a single phase slip, in 
which the number of nodes in the standing-wave pattern 
increases exactly by one. With the continuing increase 
of the control parameter g the secondary pattern eventu- 
ally undergoes another Eckhaus transition to fc2 — 6.60 
(to = 69), followed by a further Eckhaus transition to 
fca ^ 10.30 (to = 70). 

Upon decreasing the value of the control parameter g 
back to zero, the k^ pattern remains stable down to its 
saddle-node bifurcation at a value of g just below 20. As 
we further decrease g, the fc2 wave number is skipped 
and the fci wave number appears, even though the con- 
trol parameter is varied quasistatically. This transition 
is discussed in detail below. The single-mode pattern of 
wave number fci eventually reaches its saddle-node bifur- 
cation value and is replaced by the fco pattern. 

This sequence of secondary transitions, which is ex- 
pected for such a quasistatic upward sweep of the control 
parameter followed by a quasistatic downward sweep, is 
verified numerically in Fig. [21 Solid curves show the an- 
alytical values PH)) of the amplitudes b^ of the modes fc„ 
(n = 0...3), plotted in the region in which the corre- 
sponding single-mode solutions ([7]) are stable. Superim- 
posed symbols show the numerical solution of both the 
original equations of motion ^ for iV = 92 resonators. 
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FIG. 2: (Color online) The amplitude of single-mode oscil- 
lations as a function of the reduced drive amplitude g. The 
parameters are the same as in Fig. [T] Solid lines show the 
analytical values (|10p of the amplitudes h\ of the modes fc„ 
(for n = 0, . . . , 3). Note that the fco mode bifurcates super- 
critically, whereas all the other modes start at a saddle-node 
bifurcations. All modes terminate at the values of g for which 
they become Eckhaus unstable. Symbols show numerical cal- 
culations, where blue is used for upward sweeps of g and red 
is used for downward sweeps, as follows: (a) -|-s and Ds show 
upward and downward sweeps, respectively, of the original 
equations of motion ([1} of the coupled resonators, (b) As 
and ys show upward and downward sweeps, respectively, of 
the BCL amplitude equation ((6]|. (c) *s and os show upward 
and downward sweeps, respectively, of the truncated mode 
expansion equations lfT6]l for the seven modes 6_3 to 63. 



and the BCL amplitude equation ([6|), for a quasistatic 
sweep of the control parameter from 5 = up to g = 85 
and back down to g = Q. We note that in order to satisfy 
the boundary conditions when integrating the BCL am- 
phtude equation, both the 5^/^ and S'^/'^ terms of Eq. (H]), 
when substituted into the expression A+e"*'""^ -I- Ale**f 
in Eq. must be set to zero separately at the bound- 
aries. This yields a pair of conditions on B and its deriva- 
tive, at the boundaries, of the form 



dB 



-IQpX 



dB 



+ 1- 



(13) 
(14) 



To study the actual process of an Eckhaus transition 
as it takes place, we expand the general solution of the 
BCL amplitude equation in the linear modes of the array 



i?(e,r)=^fe„(r)e 



(15) 



where kn is defined in ([9]), as was done, for example, 
in a similar situation by Kramer et al. [2J|. Substitut- 
ing a truncated mode expansion psp containing a finite 
number of modes around into the BCL amplitude 
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FIG. 3: (Color online) Time evolution of the amplitudes of 
the four largest modes that participate in the Eckhaus tran- 
sition from the initial fco pattern to the k\ pattern, obtained 
by a numerical integration of the seven truncated mode equa- 
tions (|16p . for modes 6-3 to 63, using the same parameters as 
in Fig. [2] The value of the control parameter is changed from 
(? = 10 to g = 11 at r = 0, causing the initial fco pattern to be- 
come unstable. The decay of the amplitude bo is followed by 
the rise of 61 to its expected steady-state value (|10|l , where it is 
clearly seen that during the transition other modes — including 
the unstable mode — have a non-zero amplitude. 



equation ([6]), allows us to replace this partial differen- 
tial equation with a finite number of ordinary differential 
equations for the coupled mode amplitudes, 



^ = {g~kl)hn 



2j2(k,-l- 



P^m+p—n 



m,l^p,r 



(16) 



To satisfy the boundary conditions, as mentioned above 
for the single- mode solution ([7]), we take each mode am- 
plitude to be zero at the boundaries by setting all the 
phases in Eq. (fT5|) to 7r/4, and take the amplitudes 
bn to be real, keeping in mind that they can be either 
positive or negative. Note that if all mode amplitudes 
except &o are set to zero we obtain a single equation with 
n — m=p — I — r — 0, whose steady-state solution is 
the same as the single- mode solution of BCL (|10|) . 

As shown in Fig. [21 we can capture the sequence of Eck- 
haus transitions from the fcp pattern up to the ^3 pattern, 
and back down to the fco pattern through the saddle- 
nodes, by integrating seven coupled ordinary differential 
equations (calculated symbolically using MATLAB) 
that correspond to a truncated mode expansion con- 
taining the seven modes from n = —3 to n = 3. The re- 
sults agree very well with those obtained by integrating 
the full BCL amplitude equation as well as the original 
equations of motions for the resonators, yet the solution 




FIG. 4: (Color online) Time evolution of the amplitudes of the 
four largest modes as the control parameter is changed from 
g = 20 to g = 19, obtained by a numerical integration of the 
seven truncated mode equations (|16|) . for the amplitudes fo_3 
to foa, using the same parameters as in Fig.[T] As the value of g 
drops below the saddle node value of the ks wave number, its 
amplitude drops abruptly to zero. Then, the smallest possible 
wave number which has the largest linear growth rate over the 
zero solution, fco, grows to reach its steady-state value (I10|l . 
Nevertheless, after a short transient the fco pattern decays 
through an Eckhaus instability and the fci pattern grows to 
its steady-state value. 

in terms of a truncated mode expansion allows us to in- 
spect the transitions between patterns in greater detail. 

We take a closer look at the transient behavior during 
the first Eckhaus transition from the initial fco pattern to 
the fci pattern by plotting the time evolution of the four 
largest modes during this Eckhaus transition, as shown in 
Fig. [3l One can observe the decay of the unstable mode 
amplitude bo followed by the growth of bi to its steady- 
state value. One can also see that during the transient 
the amplitude of the unstable mode 6_i becomes non- 
zero. Its participation in the Eckhaus transition from 
the fco pattern to the fci pattern is essential, as can be 
verified by considering these two modes alone in a trun- 
cated expansion. Limiting the expansion to bo and &i 
suppresses the Eckhaus transition, and the fco pattern 
remains stable as g exceeds its expected value for the 
Eckhaus instability. The Eckhaus transition is observed 
only when the fc_i mode is included as well, correspond- 
ing to the stability calculation, performed earlier for the 
state given by Eq. (fTT|) . 

One might naively expect that the same mechanism 
causes the transition from the fca pattern to the fci pat- 
tern at g = 19 through a double phase slip, however, 
this is not the case. Fig. 2] reveals the transient pro- 
cesses on a downward sweep of g just below the saddle 
node at 5 = 19. As g crosses the saddle node bifurca- 
tion point, the amplitude 63 drops abruptly to zero. As 
can be seen from Eq. (|16p , in the zero displacement state 
the linear growth rates of the solutions (fT5)) are g — k^, 
so the fco pattern has the largest possible growth rate 



FIG. 5: (Color online) The linear growth rate ^g,ko{Q) a-s 
a function of the wave number shift Q, plotted for different 
values of the control parameter g. The horizontal axis labels 
Q in units of AQjv, which for a finite system of A'^ = 500 
resonators, with fco — —0.075, has the value of AQn — 0.69 
(all other parameters are the same as in Fig.[l|. The Eckhaus 
instability of the initial fco pattern occurs for these parameters 
at y ~ 5.13. 



and it out-grows the other modes until its amplitude ap- 
proaches the steady state value ([TU)) . However, according 
to the eigenvalue ^2]) , at this value of g the fco pattern is 
Eckhaus unstable with respect to the fci pattern — notice 
the characteristic evolution of the modes around r = 3 
in Fig. 21 corresponding to the Eckhaus instability (cf. 
around r = 25 in Fig. [3]). Thus the fci mode is ultimately 
the selected pattern. 



V. ABRUPT CHANGE OF THE CONTROL 
PARAMETER 

For a quasistatic increase of g the Eckhaus instability 
leads to a single phase slip event and a jump by one of the 
mode number. This is no longer always the case for more 
rapid variations of g [29] . In this section we consider the 
question of pattern selection after an abrupt jump in g, so 
that single-mode states of different wave numbers com- 
pete with each other after the system is initiated in an 
Eckhaus unstable state. We consider a scenario in which 
the system is initiated in the fco single-mode state (O , af- 
ter which the control parameter g is abruptly increased so 
that the fco wave number is no longer stable, while many 
other wave numbers become simultaneously stable. In 
order to predict the single-mode pattern that is selected 
we use our previous expression (|12p for the eigenvalue 
^g.ko (Q) to calculate the linear growth rate of perturba- 
tions of single-mode patterns of wave number ko + Q. In 
Fig. [5] we plot the Xg^kg (Q) as a function of Q for four 
different values of g, illustrating the dependence of the 
fastest growing wave number fco + Qmax on g. The wave 
number fco + Qmax with the largest linear growth rate 
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predictions for the selected wave numbers are verified nu- 
merically by integrating the BCL amplitude equation 
as well as the original equations of motion ([T]) of the cou- 
pled resonators. We note that for the parameters used 
the stability balloon contains about 10 modes for each of 
the values taken for the control parameter. For g = 21, 
for example, all modes with wave numbers from /ca to 
fci6 are stable. Finally, by following the amplitude of the 
growing mode as a function of time in the numerical solu- 
tion of the BCL amplitude equation ([S]) , it is possible to 
extract the linear growth rate of the mode numerically. 
We have done so and found that the numerically calcu- 
lated growth rates agree to within 2% with the analytical 
values of Xg^ko {Q max J • 



FIG. 6: The wave number shift Qmax with the maximal 
growth rate (in units of AQjv) as a function of g. The blue 
solid line shows Qmax for the same parameters used in Fig. O 
and the green dashed line gives Qmax for an infinite sys- 
tem (|17[) . The open red circles are the wave number shifts 
that are observed numerically by integrating the BCL am- 
plitude equation The full black circles are wave number 
shifts that are obtained by a numerical integration of the orig- 
inal equations of motion ([T}, calculated for select values of g. 
The numerical calculations are initialized with the ko solution 
and random small-amplitude noise, to initiate the growth of 
competing patterns. We emphasize that the results of the nu- 
merical solution of the BCL amplitude equation (|SJ are not 
sensitive to the noise amplitude as long as it is sufficiently 
small. 

Xmax, which is expected to overcome all other modes, is 
obtained by finding the maximum of Xg^ko {Q) a-s a func- 
tion of Q, yielding 

Qlax = (3fc2-5&^^fco-36|,+26t) 
{3kl-nblko + 3bl+8bl) 

3{ko-bl)i3h-5bl) ' ^ '> 

and 

_ i3kl-5blko-3bl+2bl)^ 

3iko^bl)i3ko-5bl) ' ^^'^ 

where bkg is the steady-state amplitude of the unstable 
ko mode, as given by Eq. pH)) . that depends on the actual 
value of g. 

For a finite system the selected wave number is ex- 
pected to be the fc„ — defined in Eq. ^ — which has the 
largest linear growth rate. The Eckhaus instability is 
triggered by random small-amplitude noise. In our finite 
system the difference between growth rates of different 
modes is expected to be sufficiently large so that by the 
time nonlinear effects are important, the amplitude of the 
fastest growing mode far exceeds those of other destabi- 
lizing Eckhaus modes, and it will reach its steady state 
value. Fig. [5] shows Qmax for an infinite system and for 
a finite system of = 500 resonators, where the two 
curves should tend to one another as N — > oo. These 



VI. RAMPS OF THE CONTROL PARAMETER 

We finish by considering a scenario in which the con- 
trol parameter g varies smoothly with time — this is often 
called a control parameter ramp. To simplify the analy- 
sis we consider slow variation in time, where dg/dr <C 1, 
so that we can use the expressions ([7]) and (|10p , obtained 
earlier for the steady-state single-mode solutions of the 
BCL amplitude equation, with a simple replacement of 
the previously constant g by a time dependent g{T), 

B(e,r) =a„(T)e'('^-^-"«), (19) 
an{rf = (A:„ - 1) + ^{k„ - 1)2 + (g(r) - k^). (20) 

Thus, a„{T) would be the steady-state amplitude of the 
pattern with wave number kn if the drive were varied 
quasistatically to its instantaneous value at time t. For 
ramps that are not quasistatic we expect the actual am- 
plitude, which we denote (r), to lag behind its ex- 
pected value for a quasistatic ramp. This time lag phe- 
nomenon is known from experiments measuring the heat 
flow in a Rayleigh-Benard cell [26| . 

The specific scenario we examine is one in which the 
control parameter increases linearly in time from zero, 
g = ar with a ^ 1. Initially, the system is expected to 
evolve to the single-mode state (jl9p with wave number 
/co- As g increases, this solution becomes Eckhaus un- 
stable and a transition is expected to a different pattern, 
which eventually becomes Eckhaus unstable as well. It is 
the first of these Eckhaus transitions that we treat ana- 
lytically below, as well as test numerically using the BCL 
amplitude equation. In order to obtain interesting mode 
competition, even for a ^ 1, we increase the number of 
resonators to = 1230, thus increasing the number of 
stable single-mode solutions for any particular value of 
g. For such a number of resonators we no longer perform 
numerical calculations on the original coupled equations 
of motion ([1]). The stability balloon for N — 1230 res- 
onators is shown in Fig. [T] Due to the large number of 
resonators, the Eckhaus boundaries for infinite and fi- 
nite systems are almost the same. The dashed vertical 
lines mark the values of possible wave numbers fc„ for 
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FIG. 7: (Color online) Stability balloon for an array of 
A'^ = 1230 resonators and the same parameters as in Fig. [T] 
which yields kg ~ 0.075 and AQjv — 0.28. The values of k„ 
for n = — 1...10 are marked with vertical dotted lines. The 
thick dashed line is the neutral stability curve g = k^, and 
the segment in which it is the lower boundary \k\ < A(5jv/2 
is marked by a solid red line. The black dotted and solid 
curves (which are almost indistinguishable) are the Eckhaus 
boundaries for an infinite and a finite system respectively. 
The solid green line is the saddle node g — 2k ~ 1 which is 
the lower boundary of oscillations for k > 2.47. 
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FIG. 8: (Color online) The three relevant Fourier amplitudes 
of the numerical solution of the amplitude equation (|6} for g = 
10~''r and the same parameters as in Fig. [T] The quasistatic 
values of the amplitudes (|20|l are plotted in thin black lines. 
For a — 10~* we expect a double phase slip from the kg 
pattern to the ^2 pattern as can be inferred from Fig. [TT] 
The inset demonstrates the time lag at early times between 
the actual amplitude of the ko mode and its quasistatic value. 



n = — 1...10. For ri = (fc ~ 0.075) the lower boundary 
of oscillations is the neutral stability curve, for n — 1...8 
the lower boundary is the Eckhaus instability curve, and 
for n > 9 (fc ~ 2.5) it is the saddle node line. 

A typical response of the system to a linear ramp of 
the drive amplitude is shown in Fig. [5] for a ramp rate 



FIG. 9: (Color online) The real part of B(^, r), obtained by a 
numerical integration of the BCL amplitude equation ([6]) for 
g — 10~^r, using the same parameters as in Fig. [7] The initial 
zero-state 5(^,0) = evolves into the ko state, which then 
undergoes a sequence of Eckhaus transitions as g increases 
in time — the first two transitions involve single phase slips, 
while the third involves a double phase slip. 



of a = lO^"'. One clearly sees the amplitude of the ko 
mode, which forms initially from the zero-state, becom- 
ing Eckhaus unstable around g — 7 and undergoing a 
double phase slip to the k2 mode. Thin black lines show 
the quasistatic values ([20| of the amplitudes of these two 
modes as a function of ^(t), while the blue dot-dashed 
curve and the green dashed curve show the actual values 
of these two amplitudes as obtained by Fourier trans- 
forming the numerical solution of the BCL equation. The 
curves are distinguishable from each other only at very 
early times, shown in the inset of Fig. [51 clearly demon- 
strating the time it takes the actual amplitude ao(r) to 
"catch up" with the quasistatic value ao ('''), from the 
zero displacement state. After this initial time lag the 
system responds sufficiently quickly so that the ramp be- 
comes effectively quasistatic. The only points where the 
time dependence of the ramp is still evident are the Eck- 
haus instability points at which different ramp rates are 
expected to lead to different transitions. A typical se- 
quence of such transitions is shown in Fig. [5] for a ramp 
rate of a = 10^^. 

To analytically predict the first Eckhaus transition, 
given the ramp rate a, it is useful to introduce a more 
compact notation for the eigenvalues (|12p. which now de- 
pend on time [30^, denoting Ag(T-) (nAQAr) = Xnlgij)]. 
The first five of these eigenvalues are plotted in Fig. [TO] 
as a function of ^(r). Note that as n increases, the corre- 
sponding eigenvalue A„ becomes positive at a later point 
in time, which we denote as r„, but grows more rapidly 
than the smaller-n eigenvalues. At time r„ the amplitude 
of the kn mode is expected to start growing from its ini- 
tial value a„(r„), which in a real physical system is set by 
the noise floor. In our analysis below we take this initial 
value to be the same as the accuracy of the numerical rou- 
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FIG. 10: The first five eigenvalues A,i, plotted as a function of 
g = ar using the parameters of Fig. [7] As Ai turns positive 
the ko solution becomes Eckhaus unstable, but the selected 
pattern depends on the ramp rate a as explained in the text. 



tine that is used for integrating the BCL equation. Once 
the kn pattern starts growing it competes with all the 
other single- mode patterns with positive growth rates. 
We expect the pattern that is eventually selected to be 
the one whose amplitude is first to reach the quasistatic 
value a„(r), given by Eq. ((20|) . Thus, lower- ti modes 
have an advantage for small ramp rates a because their 
growth rates become positive earlier. At higher ramp 
rates, due to the time-lag phenomenon shown above, the 
higher-n modes have an advantage because their eigen- 
values increase more rapidly in time. This gives rise to 
an interesting competition between the possible stable 
patterns. A similar situation was observed in a system 
described by the stochastic time-dependent Ginzburg- 
Landau equation 

Owing to the slow ramp rates, and the fact that the 
second eigenvalue associated with each mode remains 
negative, we can estimate the growth of the n^^ ampli- 
tude, in the linear regime, from its initial value at r„ to 
be 



a„(r) = a„(r„)e' 



CT„(r,r„) 



where 



o-n(T,r„) = J Xn[g{T')]dT' 



(21) 



(22) 



A comparison of these expressions for a„(r) for the dif- 
ferent patterns allows us to determine which is the first 
to reach its quasistatic value a„(r), and provides a simple 
scheme for predicting the selected pattern following the 
Eckhaus instability of the initial ko pattern. These ana- 
lytical predictions for the selected pattern kn are shown 
as a function of the ramp rate a in Fig. Ill[ and are nicely 
verified by numerical integration of the BCL amplitude 
equation As expected, for small values of a there is 
a single phase slip to the pattern with wave number fci. 



FIG. 11: (Color online) The number of phase slips Q/AQjv, 
that are observed following the Eckhaus instability of the ini- 
tial ko pattern, plotted as a function of the ramp rate a for 
a linear ramp of the drive g = ar. Parameters are the same 
as in Fig. [T] Red circles are the actual values observed in 
the numerical integration of the BCL amplitude equation © . 
The blue line shows the predicted values from the linear anal- 
ysis described in the text, where the initial amplitude of each 
mode, when its eigenvalue becomes positive, is taken to be 
an{T„) = 10~^^, which is the accuracy of the time integration 
in the numerical solution of 



As a is further increased this changes to a double phase 
slip to the k2 pattern (as demonstrated earlier in Fig. [8]), 
followed by a transition to the ks pattern, and so on. 



VII. CONCLUSIONS 

We have investigated the sequence of single- 
mode standing-wave patterns to be expected in one- 
dimensional arrays of parametrically-driven oscillators 
for time varying drive strengths. An amplitude equa- 
tion approach provides a general treatment in terms of a 
universal stability diagram on a plot with scaled versions 
of the driving strength and wave number as axes. This 
immediately shows the type of instability that will be en- 
countered on varying parameters, and gives qualitative 
insights on the mode jumps to be expected. For exam- 
ple, for quasistatic parameter variations, we find that the 
jump in the mode number is always unity if the control 
parameter is increased so that the Eckhaus instability 
operates, but larger jumps are often seen if the control 
parameter is decreased so that a saddle node bifurcation 
occurs. For more rapid increases in the control parame- 
ter larger jumps in the mode number may also occur, and 
can be predicted simply from the eigenvalue equation for 
the Eckhaus instability. We give explicit results for exam- 
ples of an abrupt jump and a slow temporal ramp in the 
control parameter. In all cases we checked, simulations 
of the original oscillator equations of motion confirm the 
results based on the amplitude equation. 

In the Buks and Roukes experiments on 
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parametrically-driven oscillator arrays [l7| which 
motivated this study the frequency of the drive was 
swept, rather than the strength of the driving. Since the 
drive frequency is involved in setting the wave number 
of the resonant mode that goes unstable, and these two 
parameters are involved in a complicated way in the ex- 
pressions for the scaled drive and wave number variables 
(see BCL), it is more difficult to display the variation 
on the scaled stability plot of Fig. [TJ For quasistatic or 
abrupt variations this is immaterial, since the behavior is 
determined by where the stability boundary is crossed in 
the former case, or the relationship of the final point on 
the stability plot to the stability boundaries in the latter 
case. These can be estimated quite easily, so that the 
behavior for quasistatic or abrupt jumps can be readily 
predicted. In particular we expect single mode jumps for 
quasistatic variations of the drive frequency in the sense 
leading to a crossing of the Eckhaus boundary, and the 
possibility of larger mode jumps for quasistatic sweeps 
in the reverse direction. This is qualitatively consistent 
with the numerical results of Lifshitz and Cross [l^ 
for a numerical model of 67 parametrically driven 
resonators, who found more jumps in the solution on 
decreasing the frequency than on increasing it. (In the 
experiments of Buks and Roukes only upward frequency 
sweeps were performed.) A more detailed comparison 
with experiment would require a better knowledge of 
the parameters of the MEMS or NEMS devices so 
that the frequency variation could be mapped onto the 
reduced stability diagram of the amplitude equation. In 
future experiments, upward and downward sweeps of 
the strength of the driving would provide a more direct 
comparison with the theory we have developed. 
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APPENDIX A: EXPLICIT FORM OF 
SINGLE-MODE SOLUTIONS 

When substituting the single mode solution ([7]) back 
into ^ and (jH), and then into one obtains extended 
single-mode standing- wave parametric oscillations at half 
the drive frequency, whose explicit form is given by 

u{x,t) ~ e^/^(5^/''4S'f,6fe sin(gf„ia;) [cos(7r/4 - Wpt) 
-I- tan(Q!) sin(7r/4 — Wpi)] 

X cos{7t/4 ~ ojpt ~ a), 
where we have defined 



(Al) 



Unic^)^6'/^^co,sin^{^^){bl-k), (A2) 
and the scale factor 



3V3 ^ 



ujr, sm 



(f) 



(A3) 



To satisfy the boundary conditions u(0, t) = u{N+l, t) 
0, the wave numbers qm satisfy the equation 



iV+ 1 



9p 



where 



AO 



N 



1 3A2 sin(gp) TT 



(A4) 



(A5) 
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Note that for the positive square-root branch (llOp b\{k — 
1 — h\) < 0, implying that Ag,fc(0) — 0, while for the 



negative square-root solution 



bl{k 



I - hi) > 0, im- 



plying that Ag_fc(0) > 0. As a consequence the positive 
square-root solution (|10p is stable with respect to small 
perturbations with the same wave number k, while the 
negative square-root solution is unstable. 
The behavior when the saddle node instability is encoun- 
tered on decreasing g is expected to be similar for slow 
or rapid variation of g, since the dynamics is simply the 
decay of the mode to zero, with the subsequent growth 
of other modes. 

For slow ramps the instantaneous rate of growth of per- 
turbations is well approximated by ignoring the time de- 
pendence of the parameter. 



